
function [f,p,nk,namevec,css,nss,cy,xss,yss,A,B]=IBCmodel(params)



jj=1;
chi                 = params(jj); jj=jj+1;
xi                  = params(jj); jj=jj+1;
xiw                  = params(jj); jj=jj+1;
rhog                = params(jj); jj=jj+1;
sigma               = params(jj); jj=jj+1;
gamma               = params(jj); jj=jj+1;
marketscomplete     = params(jj); jj=jj+1;
adjcostCEE          = params(jj); jj=jj+1;
fixrate             = params(jj); jj=jj+1;
psi_t                = params(jj); jj=jj+1;
psi_g                = params(jj); jj=jj+1;
rhoi                = params(jj); jj=jj+1;
phipi               = params(jj); jj=jj+1;
phitg               = params(jj); jj=jj+1;
psiy                = params(jj); jj=jj+1;
gy                  = params(jj);jj=jj+1;
n                   = params(jj);jj=jj+1;
hab                 = params(jj);jj=jj+1;
phigap               = params(jj);jj=jj+1;
rho_comp             = params(jj);jj=jj+1;
cgg00             = params(jj);jj=jj+1;
noinv             = params(jj);jj=jj+1;
nolongrate             = params(jj);jj=jj+1;
importshare         = params(jj);jj=jj+1;
distax             = params(jj);jj=jj+1;
cpitarg  = params(jj);jj=jj+1;
avtarg  = params(jj);

%% Other parameters
mu          = 3 ;           % inverse Frisch

if distax
    tau         = gy;            % tax rate in steady state
elseif ~distax
    tau = 0;
    
end
    
rhoz        = .6;           % persistence of technology shocks

if n < 1
    omega = importshare/(1-n)/(1-gy);
elseif n == 1
    omega = 1;    
end    


%% STEADY STATE TARGETS
nss         = .33;            % hours in steady state
epsilon     = 11;             % elasticity of substitution between intermediate goods
nu          = 11;
ky          = 10;             % capital-output ratio (quarterly)
xy          = 0.24;           % investment ratio (including durable consumption)
labshare    = 0.55;            
delta       = xy/ky;          % depreciation rate
theta       = 1-epsilon/(epsilon-1)*labshare;  % labor share
beta        = 1/((epsilon-1)/epsilon*(1-tau)*theta/ky+1-delta); % euler equation in steady state
cy          = 1-delta*ky-gy;
yss         = ky^(theta/(1-theta))*nss;
css         = yss*cy;
xss         = yss*delta*ky;

% disp(['delta: ' num2str(xy)])
% disp(['cy: ' num2str(cy)])
% disp(['beta: ' num2str(beta)])
% disp(['xy: ' num2str(xy)])
% disp(['theta: ' num2str(theta)])
vartheta = (1-tau)*((nu-1)/nu)*((epsilon-1)/epsilon)*(1-theta)*yss*css^(-gamma)/(nss^(1+mu));

% disp(['vartheta: ' num2str(vartheta)])

%% INITIALIZE VARIABLES
jkl =1;
ga1        = jkl; namevec=['ga1'];                      jkl = jkl +1;
ga2        = jkl; namevec=strvcat(namevec,'ga2');       jkl = jkl +1;
z1         = jkl; namevec=strvcat(namevec,'z1');        jkl = jkl +1;
z2         = jkl; namevec=strvcat(namevec,'z2');        jkl = jkl +1;
nu1        = jkl; namevec=strvcat(namevec,'nu1');       jkl = jkl +1;
nu2        = jkl; namevec=strvcat(namevec,'nu2');       jkl = jkl +1;
k1         = jkl; namevec=strvcat(namevec,'k1');        jkl = jkl +1;
k2         = jkl; namevec=strvcat(namevec,'k2');        jkl = jkl +1;
x1lag      = jkl; namevec=strvcat(namevec,'x1lag');                        jkl = jkl +1;
x2lag      = jkl; namevec=strvcat(namevec,'x2lag');                        jkl = jkl +1;
i1lag      = jkl; namevec=strvcat(namevec,'i1lag');                        jkl = jkl +1;
i2lag      = jkl; namevec=strvcat(namevec,'i2lag');                        jkl = jkl +1;
totlag     = jkl; namevec=strvcat(namevec,'totlag');                        jkl = jkl +1;
fbond1        = jkl; namevec=strvcat(namevec,'Net foreign assets');                        jkl = jkl +1;
d1         = jkl; namevec=strvcat(namevec,'Public debt');                        jkl = jkl +1;
d2         = jkl; namevec=strvcat(namevec,'Public debt^*');                        jkl = jkl +1;
qb1lag     = jkl; namevec=strvcat(namevec,'qb1lag');                        jkl = jkl +1;
qb2lag     = jkl; namevec=strvcat(namevec,'qb2lag');                        jkl = jkl +1;
y1lag     = jkl; namevec=strvcat(namevec,'y1lag');                        jkl = jkl +1;
y2lag     = jkl; namevec=strvcat(namevec,'y2lag');                        jkl = jkl +1;
c1lag     = jkl; namevec=strvcat(namevec,'c1lag');                        jkl = jkl +1;
c2lag     = jkl; namevec=strvcat(namevec,'c2lag');                        jkl = jkl +1;
w1lag     = jkl; namevec=strvcat(namevec,'w1lag');                        jkl = jkl +1;
w2lag     = jkl; namevec=strvcat(namevec,'w2lag');                        jkl = jkl +1;
tau1       = jkl; namevec=strvcat(namevec,'Tax rate');                        jkl = jkl +1;
tau2       = jkl; namevec=strvcat(namevec,'tau2');                        jkl = jkl +1;
qa1lag       = jkl; namevec=strvcat(namevec,'qa1lag');                        jkl = jkl +1;
g1lag      = jkl; namevec=strvcat(namevec,'g1lag');                        jkl = jkl +1;
comp1lag       = jkl; namevec=strvcat(namevec,'comp1lag');                        jkl = jkl +1;
comp2lag      = jkl; namevec=strvcat(namevec,'comp2lag');                        jkl = jkl +1;
e0      = jkl; namevec=strvcat(namevec,'e0');                        jkl = jkl +1;
g2lag      = jkl; namevec=strvcat(namevec,'g2lag');     

nk = jkl;   %number of predetermined variables

jkl = jkl +1;

gtot1        = jkl; namevec=strvcat(namevec,'Total spending');                        jkl = jkl +1;
gtot2        = jkl; namevec=strvcat(namevec,'Total spending^*');                        jkl = jkl +1;

comp1        = jkl; namevec=strvcat(namevec,'comp1');                        jkl = jkl +1;
comp2        = jkl; namevec=strvcat(namevec,'comp2');                        jkl = jkl +1;
n1        = jkl; namevec=strvcat(namevec,'n1');                        jkl = jkl +1;
n2        = jkl; namevec=strvcat(namevec,'n2');                        jkl = jkl +1;
ds        = jkl; namevec=strvcat(namevec,'ds');                        jkl = jkl +1;
tot       = jkl; namevec=strvcat(namevec,'Terms_of_trade');                        jkl = jkl +1;
y1        = jkl; namevec=strvcat(namevec,'Output');                        jkl = jkl +1;
y2         = jkl; namevec=strvcat(namevec,'Output^*');                        jkl = jkl +1;
qk1        = jkl; namevec=strvcat(namevec,'qk1');                        jkl = jkl +1;
qk2        = jkl; namevec=strvcat(namevec,'qk2');                        jkl = jkl +1;
tbpos         = jkl; namevec=strvcat(namevec,'Trade balance');                        jkl = jkl +1;
x1          = jkl; namevec=strvcat(namevec,'x1');                        jkl = jkl +1;
x2         = jkl; namevec=strvcat(namevec,'x2');                        jkl = jkl +1;
i1         = jkl; namevec=strvcat(namevec,'Policy rate');                        jkl = jkl +1;
i2         = jkl; namevec=strvcat(namevec,'i2');                        jkl = jkl +1;
pia1       = jkl; namevec=strvcat(namevec,'Inflation');                        jkl = jkl +1;
pib2       = jkl; namevec=strvcat(namevec,'pib2');                        jkl = jkl +1;
mcr1       = jkl; namevec=strvcat(namevec,'mcr1');                        jkl = jkl +1;
rx         = jkl; namevec=strvcat(namevec,'Real exchange rate');                        jkl = jkl +1;
t1         = jkl; namevec=strvcat(namevec,'Taxes');                        jkl = jkl +1;
t2         = jkl; namevec=strvcat(namevec,'Taxes^*');                        jkl = jkl +1;
r1         = jkl; namevec=strvcat(namevec,'r1');                        jkl = jkl +1;
r2         = jkl; namevec=strvcat(namevec,'r2');                        jkl = jkl +1;
w1         = jkl; namevec=strvcat(namevec,'Real wage');                        jkl = jkl +1;
w2         = jkl; namevec=strvcat(namevec,'w2');                        jkl = jkl +1;
mcr2       = jkl; namevec=strvcat(namevec,'mcr2');                        jkl = jkl +1;
lam1       = jkl; namevec=strvcat(namevec,'lam1');                        jkl = jkl +1;
lam2       = jkl; namevec=strvcat(namevec,'lam2');                        jkl = jkl +1;
beta1      = jkl; namevec=strvcat(namevec,'beta1');                        jkl = jkl +1;
beta2      = jkl; namevec=strvcat(namevec,'beat2');                        jkl = jkl +1;
pi1        = jkl; namevec=strvcat(namevec,'CPI inflation');                        jkl = jkl +1;
pi2        = jkl; namevec=strvcat(namevec,'pi2');                        jkl = jkl +1;
test       = jkl; namevec=strvcat(namevec,'test');                        jkl = jkl +1;
c1         = jkl; namevec=strvcat(namevec,'c1');                        jkl = jkl +1;
c2         = jkl; namevec=strvcat(namevec,'c2');                        jkl = jkl +1;
g1         = jkl; namevec=strvcat(namevec,'g1');                        jkl = jkl +1;
g2         = jkl; namevec=strvcat(namevec,'g2');                        jkl = jkl +1;
qb1        = jkl; namevec=strvcat(namevec,'qb1');                        jkl = jkl +1;
qa2        = jkl; namevec=strvcat(namevec,'qa2');                        jkl = jkl +1;
qa1        = jkl; namevec=strvcat(namevec,'qa1');                        jkl = jkl +1;
qb2        = jkl; namevec=strvcat(namevec,'qb2');                        jkl = jkl +1;
rr1         = jkl; namevec=strvcat(namevec,'Real interest rate');                        jkl = jkl +1;
rr2         = jkl; namevec=strvcat(namevec,'rr2');                        jkl = jkl +1;
x_y         = jkl; namevec=strvcat(namevec,'Investment');                        jkl = jkl +1;
c_y         = jkl; namevec=strvcat(namevec,'Consumption');                        jkl = jkl +1;
g_y        = jkl; namevec=strvcat(namevec,'Government spending');                        jkl = jkl +1;
c_y2       = jkl; namevec=strvcat(namevec,'Consumption2');                        jkl = jkl +1;
x_y2       = jkl; namevec=strvcat(namevec,'Investment ratio EMU');                        jkl = jkl +1;
piw1       = jkl; namevec=strvcat(namevec,'wage inflation');                        jkl = jkl +1;
piw2      = jkl; namevec=strvcat(namevec,'wage inflation^*');                        jkl = jkl +1;
pdv      = jkl; namevec=strvcat(namevec,'PDV');                        jkl = jkl +1;
rrl      = jkl; namevec=strvcat(namevec,'Long-term real rate');                        jkl = jkl +1;
avpi1     = jkl; namevec=strvcat(namevec,'Average inflation');                        jkl = jkl +1;     
avpi2     = jkl; namevec=strvcat(namevec,'Average inflation');                        jkl = jkl +1;     
k_y        = jkl; namevec=strvcat(namevec,'ky');                         


A = zeros(jkl,jkl);  B = zeros(jkl,jkl);  


%% household FOC

% lagrange multiplier
eqnum = 1;
B(eqnum,lam1)  = 1-hab;
B(eqnum,c1)    = gamma;
B(eqnum,c1lag)    = -hab*gamma;

eqnum = eqnum + 1;
B(eqnum,lam2)  = 1;
B(eqnum,c2)    = gamma;
B(eqnum,c2lag)    = -hab*gamma;

% labor supply
eqnum = eqnum +1;
B(eqnum,piw1) = xiw*(1+mu*nu);
A(eqnum,piw1) = beta*xiw*(1+mu*nu);
B(eqnum,lam1) = (1-xiw)*(1-xiw*beta);
B(eqnum,n1)   = -mu*(1-xiw)*(1-xiw*beta);
B(eqnum,tau1) = -tau/(1-tau)*(1-xiw)*(1-xiw*beta);
B(eqnum,w1)   = (1-xiw)*(1-xiw*beta);

eqnum = eqnum +1;
B(eqnum,piw2) = xiw*(1+mu*nu);
A(eqnum,piw2) = beta*xiw*(1+mu*nu);
B(eqnum,lam2)   = (1-xiw)*(1-xiw*beta);
B(eqnum,n2)   =  -mu*(1-xiw)*(1-xiw*beta);
B(eqnum,tau2) = - tau/(1-tau)*(1-xiw)*(1-xiw*beta);
B(eqnum,w2)   = (1-xiw)*(1-xiw*beta);

% wage inflation
eqnum = eqnum +1;
B(eqnum,piw1) = -1;
B(eqnum,w1)  = 1;
B(eqnum,w1lag) = -1;
B(eqnum,pi1)   = 1;

eqnum = eqnum +1;
B(eqnum,piw2) = -1;
B(eqnum,w2)  = 1;
B(eqnum,w2lag) = -1;
B(eqnum,pi2)   = 1;

% Euler equations
eqnum = eqnum + 1;
B(eqnum,lam1)   = 1;
B(eqnum,beta1)  = -1;
B(eqnum,i1)     = -1;
A(eqnum,lam1)   = 1;
A(eqnum,pi1)    = -1;

eqnum = eqnum+1;
B(eqnum,lam2)   = 1;
B(eqnum,beta2)  = -1;
B(eqnum,i2)     = -1;
A(eqnum,lam2)   = 1;
A(eqnum,pi2)    = -1;

if adjcostCEE
    
    % foc investment
    
    if noinv
        eqnum = eqnum+1;
        B(eqnum, x1)    = 1;
      
        eqnum = eqnum+1;
        B(eqnum, x2)    = 1;
      
    elseif ~noinv
        % foc investment
        eqnum = eqnum+1;
        B(eqnum, qk1)    =  1;
        B(eqnum, lam1)  =  -1;
        B(eqnum, x1)    = -(1+beta)*chi;
        B(eqnum, x1lag) = chi;
        A(eqnum, x1)    = -beta*chi;

        eqnum = eqnum+1;
        B(eqnum, qk2)    =  1;
        B(eqnum, lam2)  =  -1;
        B(eqnum, x2)    = -(1+beta)*chi;
        B(eqnum, x2lag) = chi;
        A(eqnum, x2)    = -beta*chi;
    end
    
    
    % foc capital
    eqnum = eqnum+1;
    B(eqnum, qk1)    = 1;
    B(eqnum, beta1) = -1;
    A(eqnum, lam1)  = (1-beta*(1-delta));
    A(eqnum, r1)    = (1-beta*(1-delta));
    A(eqnum, tau1)  = -(1-beta*(1-delta))*tau/(1-tau);    
    A(eqnum, qk1)    = beta*(1-delta);

    eqnum=eqnum+1;
    B(eqnum,qk2)    = 1;
    B(eqnum,beta2) = -1;
    A(eqnum,lam2)  = (1-beta*(1-delta));
    A(eqnum,r2)    = (1-beta*(1-delta));
    A(eqnum,tau2)  = -(1-beta*(1-delta))*tau/(1-tau);    
    A(eqnum,qk2)    = beta*(1-delta);

elseif ~adjcostCEE
    
    % foc investment
    eqnum=eqnum+1;
    B(eqnum,qk1)    = 1;
    B(eqnum,lam1)  = -1;
    B(eqnum,x1)    = -chi;
    B(eqnum,k1)    = chi;        

    eqnum=eqnum+1;
    B(eqnum,qk2)    = 1;
    B(eqnum,lam2)  = -1;
    B(eqnum,x2)    = -chi;
    B(eqnum,k2)    = chi;

    % foc capital
    eqnum=eqnum+1;
    B(eqnum, qk1)    = 1;
    B(eqnum, beta1) = -1;
    A(eqnum, lam1)  = (1-beta*(1-delta));
    A(eqnum, r1)    = (1-beta*(1-delta));
    A(eqnum, tau1)  = -(1-beta*(1-delta))*tau/(1-tau);    
    A(eqnum, qk1)    = beta*(1-delta);
    A(eqnum, x1)    = beta*delta*chi;
    A(eqnum, k1)    = -beta*delta*chi;

    eqnum=eqnum+1;
    B(eqnum,qk2)    = 1;
    B(eqnum,beta2) = -1;
    A(eqnum,lam2)  = (1-beta*(1-delta));
    A(eqnum,r2)    = (1-beta*(1-delta));
    A(eqnum,tau2)  = -(1-beta*(1-delta))*tau/(1-tau);    
    A(eqnum,qk2)    = beta*(1-delta);
    A(eqnum,x2)    = beta*delta*chi;
    A(eqnum,k2)    = -beta*delta*chi;
end

% law of motion capital
eqnum=eqnum+1;
A(eqnum,k1)    = 1;
B(eqnum,k1)    = 1-delta;
B(eqnum,x1)  = delta;

eqnum=eqnum+1;
A(eqnum,k2)    = 1;
B(eqnum,k2)    = 1-delta;
B(eqnum,x2)   = delta;


%% intermediate good firms

% philipps curve
eqnum=eqnum+1;
A(eqnum,pia1)  = beta*xi/(1-xi)/(1-xi*beta);
B(eqnum,pia1)   = xi/(1-xi)/(1-xi*beta); 
B(eqnum,mcr1)  = -1;
B(eqnum,qa1)   = 1;

eqnum=eqnum+1;
A(eqnum,pib2)  = beta*xi/(1-xi)/(1-xi*beta);
B(eqnum,pib2)  = xi/(1-xi)/(1-xi*beta); 
B(eqnum,mcr2)  = -1;
B(eqnum,qb2)   = 1;


% production function
eqnum=eqnum+1;
B(eqnum,y1)    = -1;
B(eqnum,z1)    = 1;
B(eqnum,n1)    = 1-theta;
B(eqnum,k1)    = theta;

eqnum=eqnum+1;
B(eqnum,y2)    = -1;
B(eqnum,z2)    = 1;
B(eqnum,n2)    = 1-theta;
B(eqnum,k2)    = theta;

% real marginal costs
eqnum=eqnum+1;
B(eqnum,mcr1)  =-1;
B(eqnum,r1)    =1;
B(eqnum,z1)    =-1;
B(eqnum,k1)    =1-theta;
B(eqnum,n1)    =-(1-theta);

eqnum=eqnum+1;
B(eqnum,mcr2)  =-1;
B(eqnum,r2)    =1;
B(eqnum,z2)    =-1;
B(eqnum,k2)    =1-theta;
B(eqnum,n2)    =-(1-theta);

% optimal factor input
eqnum=eqnum+1;
B(eqnum,r1)    = 1;
B(eqnum,k1)    = 1;
B(eqnum,w1)    = -1;
B(eqnum,n1)    = -1;

eqnum=eqnum+1;
B(eqnum,r2)    = 1;
B(eqnum,k2)    = 1;
B(eqnum,w2)    = -1;
B(eqnum,n2)    = -1;

% relative prices
eqnum=eqnum+1;
B(eqnum,rx)=1;
B(eqnum,qb1)=1;
B(eqnum,qb2)=-1;

eqnum=eqnum+1;
B(eqnum,rx)=1;
B(eqnum,qa1)=1;
B(eqnum,qa2)=-1;

% defines ds - from first difference of loop
eqnum=eqnum+1;
B(eqnum,qb1) =-1;
B(eqnum,qb1lag) =1;
B(eqnum,ds) =1;
B(eqnum,pi2) = 1;
B(eqnum,pi1) = -1;
B(eqnum,qb2) =1;
B(eqnum,qb2lag) =-1;

eqnum=eqnum+1;
B(eqnum,qa1)=1-(1-n)*omega;
B(eqnum,qb1)=(1-n)*omega;     

eqnum=eqnum+1;
B(eqnum,qb2)=1-n*omega;
B(eqnum,qa2)=n*omega;



%% monetary policy

if fixrate 
    eqnum=eqnum+1;
    B(eqnum,ds) =1;
    
    eqnum=eqnum+1;
    B(eqnum,i2)   =-1;
    B(eqnum,i2lag)= rhoi;
    B(eqnum,pi2)  = .5*phipi*(1-rhoi);
    B(eqnum,pi1)  = .5*phipi*(1-rhoi);
    B(eqnum,nu2)  = 1;
    
elseif ~fixrate
    eqnum=eqnum+1;
    B(eqnum,i1)   =-1;
    B(eqnum,i1lag)= rhoi;
   if cgg00
        A(eqnum,pia1)  = -phipi*(1-rhoi);
        A(eqnum,y1)  =   -phigap*(1-rhoi);
    elseif ~cgg00
        if ~cpitarg
            if avtarg
                B(eqnum,avpi1)  = phipi*(1-rhoi);
            elseif ~avtarg
                B(eqnum,pia1)  = phipi*(1-rhoi);
            end
        elseif cpitarg
            B(eqnum,pi1)  = phipi*(1-rhoi);
        end
        B(eqnum,y1)  =   phigap*(1-rhoi);
    end
    B(eqnum,nu1)  = 1;
    
    eqnum=eqnum+1;
    B(eqnum,i2)   =-1;
    B(eqnum,i2lag)= rhoi;
    if cgg00
        A(eqnum,pib2)  = -phipi*(1-rhoi);
        A(eqnum,y2)   = -phigap*(1-rhoi);
     elseif ~cgg00
         if ~cpitarg
             if avtarg
                B(eqnum,avpi2)  = phipi*(1-rhoi);
             elseif ~avtarg   
                B(eqnum,pib2)  = phipi*(1-rhoi);
             end
         elseif cpitarg
             B(eqnum,pi2)  = phipi*(1-rhoi);
         end         
         B(eqnum,y2)   = phigap*(1-rhoi);
     end    
    B(eqnum,nu2)  = 1;
end


%% fiscal policy

% budget constraint
eqnum=eqnum+1;
B(eqnum, t1)   = -1;
A(eqnum, d1)   = beta;
B(eqnum, d1)   = 1;
B(eqnum, y1)   = -tau;
B(eqnum, tau1) = -tau;
B(eqnum, qa1)  = -tau;
B(eqnum, gtot1)   = gy;
B(eqnum, qa1)  = gy;

eqnum=eqnum+1;
B(eqnum, t2)   = -1;
A(eqnum, d2)   = beta;
B(eqnum, d2)   = 1;
B(eqnum, y2)   = -tau;
B(eqnum, tau2) = -tau;
B(eqnum, qb2)  = -tau;
B(eqnum, gtot2)   = gy;
B(eqnum, qb2)  = gy;

% government spending (note g is measured in terms of percentage
% deviations from its own steady state value)
eqnum = eqnum+1;
B(eqnum,g1)     = -1;
B(eqnum,g1lag)  = rhog;
B(eqnum,y1lag)  = psiy/gy;
B(eqnum,d1)     = -psi_g/gy;
B(eqnum,ga1)    = 1;

eqnum=eqnum+1;
B(eqnum,g2)     = -1;
B(eqnum,g2lag)  = rhog;
B(eqnum,y2lag)  = psiy/gy;
B(eqnum,d2)     = -psi_g/gy;
B(eqnum,ga2)    = 1;


% Transfers (note: transfers are measured in output units)
eqnum=eqnum+1;
B(eqnum, t1)    = -1;
if ~distax
B(eqnum,d1)     = psi_t;
B(eqnum,g1)     = phitg*gy;
B(eqnum,qa1)    = phitg*gy;
end

eqnum=eqnum+1;
B(eqnum, t2)    = -1;  
if ~distax
B(eqnum,d2)     = psi_t;
B(eqnum,g2)     = phitg*gy;
B(eqnum,qb2)    = phitg*gy;
end

% tax rate 
eqnum=eqnum+1;
A(eqnum, tau1) = 1;
if distax
B(eqnum,d1)     = psi_t;
end

eqnum=eqnum+1;
A(eqnum, tau2) = 1;
if distax
B(eqnum,d2)     = psi_t;
end


%% market clearing

% intermediate good bundles
eqnum=eqnum+1;
B(eqnum,y1)    = -1;
B(eqnum,qa1)   = (gy-1)*sigma;
B(eqnum,c1)    = (1-(1-n)*omega)*cy;
B(eqnum,x1)    = (1-(1-n)*omega)*xy;
B(eqnum,c2)    = (1-n)*omega*cy;
B(eqnum,x2)    = (1-n)*omega*xy;
B(eqnum,rx)    = -(1-n)*omega*sigma*(cy+xy);
B(eqnum,gtot1)    = gy;
 
eqnum=eqnum+1;
B(eqnum,y2)    = -1;
B(eqnum,qb2)   = (gy-1)*sigma;
B(eqnum,c1)    = n*omega*cy;
B(eqnum,x1)    = n*omega*xy;
B(eqnum,rx)    = n*omega*sigma*(cy+xy);
B(eqnum,c2)    = (1-n*omega)*cy;
B(eqnum,x2)    = (1-n*omega)*xy;
B(eqnum,gtot2)    = gy;


% net exports - country 1
eqnum=eqnum+1;
B(eqnum,tbpos)    = - 1;
B(eqnum,y1)    = 1;
B(eqnum,qa1)    = (1-gy);
B(eqnum,c1)    = -cy;
B(eqnum,x1)   =  -xy;
B(eqnum,gtot1)   = -gy;


if ~marketscomplete
    
    % uip
    eqnum=eqnum+1;
    B(eqnum,i1)    = 1;
    B(eqnum,i2)    =-1;
    A(eqnum,ds)    = 1;
    
    % Domestic budget constraint: determines evolution of foreign
    % currency-bond holdings
    eqnum=eqnum+1;
    A(eqnum,fbond1)   = beta;
    A(eqnum,d1)    = beta;
    B(eqnum,c1)    = -cy;
    B(eqnum,x1)    = -xy;
    B(eqnum,y1)    = (1-tau);
    B(eqnum,qa1)   = (1-tau);       
    B(eqnum,tau1)  = - tau;
    B(eqnum,t1)    = -1;
    B(eqnum,fbond1)   = 1;
    B(eqnum,d1)    = 1;
       
    % Foreign budget constraint: test equation
    eqnum=eqnum+1;
    A(eqnum,d2)    = beta*(1-n);
    A(eqnum,fbond1)= -beta*n;
    B(eqnum,c2)    = -cy*(1-n);
    B(eqnum,x2)    = -xy*(1-n);
    B(eqnum,y2)    = (1-tau)*(1-n);
    B(eqnum,qb2)   = (1-tau)*(1-n);
    B(eqnum,tau2)  = - tau*(1-n);
    B(eqnum,t2)    = 1*(1-n);
    B(eqnum,d2)    = 1*(1-n);
    B(eqnum,fbond1)= -n;
    B(eqnum,test)  = 1;    

    % discount factor
    eqnum=eqnum+1;
    B(eqnum,beta1) = 1;
    B(eqnum,c1)    = (1-beta);
    %B(eqnum,n1)    = (1-mu)*nss/(1-nss)*(1-beta);

    eqnum=eqnum+1;
    B(eqnum,beta2) = 1;
    B(eqnum,c2)    = (1-beta);
    %B(eqnum,n2)    = (1-mu)*nss/(1-nss)*(1-beta);

elseif marketscomplete
 
    % risk sharing 
    eqnum=eqnum+1;
    B(eqnum,rx)    = 1;
    B(eqnum,lam1)  = -1;
    B(eqnum,lam2)  = 1;

    % bond holdings
    eqnum=eqnum+1;
    A(eqnum,fbond1)   = 1;

    eqnum=eqnum+1;
    B(eqnum,test) =1;

    % discount factor
    eqnum=eqnum+1;
    B(eqnum,beta1) = -1;

    eqnum=eqnum+1;
    B(eqnum,beta2) = -1;

end


%% shocks
% monetary policy
eqnum=eqnum+1;
A(eqnum,nu1)  = 1;

eqnum=eqnum+1;
A(eqnum,nu2)  = 1;
B(eqnum,nu2)  = 0;

% technology
eqnum=eqnum+1;
A(eqnum,z1)    = 1;
B(eqnum,z1)    = rhoz;

eqnum=eqnum+1;
A(eqnum,z2)    = 1;
B(eqnum,z2)    = rhoz;

% government spending
eqnum=eqnum+1;
A(eqnum, ga1)   = 1;

eqnum=eqnum+1;
A(eqnum, ga2)   = 1;

% compesation spending
eqnum=eqnum+1;
B(eqnum, comp1)   = -1;
B(eqnum, comp1lag) = rho_comp;
B(eqnum, e0)      = 1;

eqnum=eqnum+1;
B(eqnum, comp2)   = 1;
B(eqnum, comp2lag) = rho_comp;



%% definitions

% terms of trade
eqnum=eqnum+1;
B(eqnum,tot)   = -1;
B(eqnum,qa1)   = 1;
B(eqnum,qb1)   = -1;


% domestic inflation (needed for monetary policy rule)
eqnum=eqnum+1;
B(eqnum,pia1)    =-1;
B(eqnum,qa1)    =1;
B(eqnum,qa1lag)    =-1;
B(eqnum,pi1)    =1;

eqnum=eqnum+1;
B(eqnum,pib2)    =-1;
B(eqnum,qb2)    =1;
B(eqnum,qb2lag)    =-1;
B(eqnum,pi2)    =1;

% real interest rate
eqnum=eqnum+1;
B(eqnum, rr1) = -1;
B(eqnum, i1) = 1;
A(eqnum, pi1)  = 1;

% domestic long term real rate
eqnum = eqnum+1;
 B(eqnum,rrl)  = 1;
if ~nolongrate
 B(eqnum,rr1)  = -1;
 A(eqnum,rrl)  = 1;
end

eqnum=eqnum+1;
B(eqnum, rr2) = -1;
B(eqnum, i2) = 1;
A(eqnum, pi2)  = 1;


eqnum=eqnum+1;
B(eqnum,c1)    = 1;
A(eqnum,c1lag) = 1;

eqnum=eqnum+1;
B(eqnum,c2)    = 1;
A(eqnum,c2lag) = 1;

%% lags
eqnum=eqnum+1;
B(eqnum,tot)   = 1;
A(eqnum,totlag)= 1;

eqnum=eqnum+1;
B(eqnum,x1)    = 1;
A(eqnum,x1lag) = 1;

eqnum=eqnum+1;
B(eqnum,x2)    = 1;
A(eqnum,x2lag) = 1;

eqnum=eqnum+1;
B(eqnum,i1)    = 1;
A(eqnum,i1lag) = 1;

eqnum=eqnum+1;
B(eqnum,i2)    = 1;
A(eqnum,i2lag) = 1;

eqnum=eqnum+1;
B(eqnum,qb1)    = 1;
A(eqnum,qb1lag) = 1;

eqnum=eqnum+1;
B(eqnum,qb2)    = 1;
A(eqnum,qb2lag) = 1;

eqnum=eqnum+1;
B(eqnum,g1)    = 1;
A(eqnum,g1lag) = 1;

eqnum=eqnum+1;
B(eqnum,g2)    = 1;
A(eqnum,g2lag) = 1;

eqnum=eqnum+1;
B(eqnum,y1)    = 1;
A(eqnum,y1lag) = 1;

eqnum=eqnum+1;
B(eqnum,y2)    = 1;
A(eqnum,y2lag) = 1;

eqnum=eqnum+1;
B(eqnum,w1)    = 1;
A(eqnum,w1lag) = 1;

eqnum=eqnum+1;
B(eqnum,w2)    = 1;
A(eqnum,w2lag) = 1;


eqnum=eqnum+1;
B(eqnum,qa1)    = 1;
A(eqnum,qa1lag) = 1;

eqnum=eqnum+1;
B(eqnum,comp1)    = 1;
A(eqnum,comp1lag) = 1;

eqnum=eqnum+1;
B(eqnum,comp2)    = 1;
A(eqnum,comp2lag) = 1;

%% output units
eqnum=eqnum+1;
B(eqnum, c_y)       = -1;
B(eqnum, c1)        = cy;

eqnum=eqnum+1;
B(eqnum, x_y)       = -1;
B(eqnum, x1)        = xy;

eqnum=eqnum+1;
B(eqnum, g_y)       = -1;
B(eqnum, gtot1)        = gy;

eqnum=eqnum+1;
B(eqnum, k_y)       = -1;
B(eqnum, k1)        = ky;

eqnum=eqnum+1;
B(eqnum, c_y2)      = -1;
B(eqnum, c2)        = cy;

eqnum=eqnum+1;
B(eqnum, x_y2)      = -1;
B(eqnum, x2)        = xy;

% pdv tax burden
eqnum=eqnum+1;
B(eqnum, pdv)      = 1;
%B(eqnum, gtotpos)     = -(1-beta)/gy;
%B(eqnum,spos)      =  -(1-beta)*omega;
B(eqnum, t1)      = (1-beta)/gy;
A(eqnum, pdv)      = beta;

% total home g-spending
eqnum=eqnum+1;
B(eqnum, gtot1)   = -1;
B(eqnum, g1)      = 1;
B(eqnum, comp1)   = 1/gy;

% total home g-spending
eqnum=eqnum+1;
B(eqnum, gtot2)   = -1;
B(eqnum, g2)      = 1;
B(eqnum, comp2)   = 1/gy;

% innovation endowment
eqnum = eqnum+1;
A(eqnum, e0) = 1;


% average inflation

if xi == 0;
    eqnum=eqnum+1;
    B(eqnum, avpi1)   = 1;
    B(eqnum, pia1)   = -1;
    eqnum=eqnum+1;
    B(eqnum, avpi2)   = 1;
    B(eqnum, pib2)   = -1;
elseif xi>0
lambdap = (1-beta*xi)*(1-xi)/xi;
lambdaw =(1-beta*xiw)*(1-xiw)/xiw;

eqnum=eqnum+1;
B(eqnum, avpi1)   = -(lambdap+lambdaw);
B(eqnum, pia1)   = lambdap;
B(eqnum, piw1)   = lambdaw;

eqnum=eqnum+1;
B(eqnum, avpi2)   = -(lambdap+lambdaw);
B(eqnum, pib2)   = lambdap;
B(eqnum, piw2)   = lambdaw;
end


if nolongrate
    [f,p]=my_solab(A,B,nk,1.000000001); % adjustment because of permanent change in spending
elseif ~nolongrate
    [f,p]=my_solab(A,B,nk,1);
end